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ABSTRACT 

We investigate marginally stable nuclear burning on the surface of accreting neutron stars as an explanation 
for the mHz quasi-periodic oscillations (QPOs) observed from three low mass X-ray binaries. At the boundary 
between unstable and stable burning, the temperature dependence of the nuclear heating rate and cooling rate 
almost cancel. The result is an oscillatory mode of burning, with an oscillation period close to the geometric 
mean of the thermal and accretion timescales for the burning layer We describe a simple one-zone model 
which illustrates this basic physics, and then present detailed multizone hydrodynamical calculations of nuclear 
burning close to the stability boundary using the KEPLER code. Our models naturally explain the characteristic 
2 minute period of the mHz QPOs, and why they are seen only in a very narrow range of X-ray luminosities. 
The oscillation period is sensitive to the accreted hydrogen fraction and the surface gravity, suggesting a new 
way to probe these parameters. A major puzzle is that the accretion rate at which the oscillations appear in 
the theoretical models is an order of magnitude larger than the rate implied by the X-ray luminosity when 
the mHz QPOs are seen. We discuss the implications for our general understanding of nuclear burning on 
accreting neutron stars. One possibility is that the accreted material covers only part of the neutron star surface 
at luminosities Lx > lO-'^ erg s"'. 

Subject headings: accretion, accretion disks — X-rays:bursts — stars:neutron 



1. INTRODUCTION 

Low mass X-ray binaries, in which a neutron star or black 
hole accretes from a low mass companion, exhibit a range of 
periodic and quasi-periodic phenomena, ranging in frequency 
from very low frequency (mHz) noise to kHz quasi-periodic 
oscillations (QPOs) (see van der Klis 2004 for a recent re- 
view). This variability has mostly been associated with orbit- 
ing material in the accretion flow close to the compact object. 
In the case of a neutron star accretor, an important question 
is whether any of these phenomena originate from or are as- 
sociated with the neutron star surface. This is important for 
identifying the compact object as a neutron star or a black 
hole and offers a probe of the neutron star surface layers. 

Unstable nuclear burning on neutron star surfaces has been 
studied for many years, and is observed as Type I X-ray bursts. 
The accreted hydrogen (H) and helium (He) fuel accumulates 
on the surface of the star and undergoes a thin shell instability, 
giving rise to a 10-100 second burst of X-rays with typical en- 
ergy 10'''' ergs (for reviews, see Lewin, van Paradijs, & Taam 
1993, 1995; Strohmayer & Bildsten 2003). Not all sources 
show Type I X-ray bursts, however, and in many sources the 
bursts are not frequent enough to burn all of the accreted fuel 
(van Paradijs, Penninx, & Lewin 1988; in 't Zand et al. 2003). 
Bildsten (1993, 1995) suggested that a different mode of nu- 
clear burning, involving slowly propagating fires over the neu- 
tron star surface, operates at high accretion rates, and mani- 
fests itself in the power spectrum of the source as very low fre- 
quency noise (VLFN). He found an anti-correlation between 
bursting and VLFN, supporting this picture. 

Revnivtsev et al. (2001) discovered a new class of mHz 
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QPOs in three Atoll sources, 4U 1608-52, 4U 1636-53, and 
Aql X-1, which they proposed were from a special mode of 
nuclear burning on the neutron star surface rather than from 
the accretion flow. These mHz QPOs have frequencies in the 
range 7-9 mHz (timescales of 1.9-2.4 minutes). The asso- 
ciated flux variations are at the few per cent level, and are 
strongest at low photon energies (< 5 keV). This is in contrast 
to all other observed QPOs, whose amplitude generally in- 
creases with photon energy. Revnivtsev et al. (2001) showed 
that the centroid frequency of the mHz QPO was stable on 
year timescales in a given source, and the same to within tens 
of per cent in all three sources in which it was detected. In 
addition, the presence of the mHz QPOs is affected by Type 
I X-ray bursts: in 4U 1608-52, the mHz QPOs disappeared 
immediately following a Type I X-ray burst. In 4U 1608-52, 
a transient source whose luminosity is observed to change by 
orders of magnitude, the mHz QPO was only present within a 
narrow range of luminosity, Lx ~ 0.5-1 .5 x lO-'' erg s"' . 

The association of the mHz QPOs with a surface phe- 
nomenon was strengthened by the results of Yu & van der 
Klis (2002), who showed that the kHz QPO frequency is anti- 
correlated with the luminosity variations during the mHz os- 
cillation. This is opposite to the long term trend, which is 
that the kHz QPO frequency varies proportional to the X- 
ray luminosity, consistent with the inner edge of the accretion 
disk being pushed inwards at higher accretion rates. The anti- 
correlation observed by Yu & van der Klis (2002) during the 
mHz QPO cycle suggests that the inner edge of the disk moves 
outwards slightly as the luminosity increases during the cycle, 
perhaps consistent with enhanced radiation drag as the gas or- 
biting close to the neutron star is radiated by emission from 
the neutron star surface. 

The fact that the properties of the mHz QPOs suggest that 
they are associated with the neutron star surface led Revnivt- 
sev et al. (2001) to propose that a special mode of nuclear 
burning operates at luminosities Lx ~ 0.5-1.5 x lO"'^ erg s~'. 
The nature of the burning and the physics underlying the 
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characteristic timescale of w 2 minutes, however, were un- 
explained. Bildsten (1993) gave the characteristic timescales 
of the burning layer that might be associated with sub-hertz 
phenomena. The thermal timescale of the layer is ftheim ~ 10 
s, the time to accrete the fuel is tg^ca ~ 1000 s at the Edding- 
ton accretion rate, and the time for a nuclear burning "fire" 
to propagate around the surface is estimated to be ~ 10"^ s. 
Bildsten (1993) proposed that if several fires are propagating 
around the star at a given time, the time-signature would be 
broad band noise with frequencies in the mHz range. None of 
the timescales he identified match the w 2 minute mHz QPO 
period, however 

The luminosity at which mHz QPOs are observed {Lx ~ 
10" ergs-i) is significant because it is similar to the lumi- 
nosity at which a transition in burning behavior occurs, from 
frequent Type 1 X-ray bursting at low accretion rates to the 
disappearance of Type 1 X-ray bursts at high accretion rates. 
This transition is common to many X-ray bursters (e.g., Cor- 
nelisse et al. 2003), and is expected theoretically because at 
high accretion rates the fuel burns at a higher temperature, 
reducing the temperature-sensitivity of helium burning and 
quenching the thin shell instability. An outstanding puzzle is 
that theory predicts a transition accretion rate close to the Ed- 
dington rate (Bildsten 1998), which corresponds to luminosi- 
ties Lx ~ 10""^ erg s~', much larger than observed. Paczynski 
(1983) pointed out that near the transition from instability to 
stability, oscillations are expected because the eigenvalues of 
the system are complex. Narayan & Heyl (2003) extended 
Paczynski's one-zone analysis, calculating linear eigenmodes 
of truncated steady-state burning models. They too found 
complex eigenvalues near the stability boundary, and sug- 
gested that this might explain the mHz QPOs observed by 
Revnivtsev et al. (2001). The oscillation frequencies, how- 
ever, were an order of magnitude too small. 

In this paper, we show that the mHz QPO frequencies 
are, in fact, naturally explained as being due to marginally 
stable nuclear burning on the neutron star surface. At the 
boundary between unstable and stable burning, the temper- 
ature dependence of the nuclear heating rate and cooling rate 
almost cancel. The result is an oscillatory mode of burn- 
ing, with an oscillation period close to the geometric mean 
of the thermal and accretion timescales for the burning layer, 
(fthei-mfaccr)'''^ ~ 100 s. In §2, we describe a simple one-zone 
model which illustrates this basic physics, and then present 
detailed hydrodynamical calculations of nuclear burning close 
to the stability boundary in §3. We discuss the implications 
of our results in §4. In particular, if the mHz QPOs are due 
to marginally stable nuclear burning, the local accretion rate 
onto the star must be close to the Eddington rate, even though 
the global accretion rate inferred from the X-ray luminosity is 
ten times lower. 



2. A ONE-ZONE MODEL 

In this section, we discuss a simplified model of the burn- 
ing layers which illustrates the basic physics underlying the 
oscillations observed in our multizone numerical simulations. 
Following Paczynski (1983), we consider a one-zone model 
of the burning layer. The temperature, T, and thickness of the 
fuel layer, y (which we measure as a column depth, in units of 
mass per unit area), obey 
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(equivalent to eqns. [8] of Paczynski 1983). Equation ([!) de- 
scribes the heat balance, including heating of the layer by nu- 
clear reactions e, and radiative cooling -V- F/p = dF /dy ~ 
F jy. The heat capacity at constant pressure is cp, and F is 
the outwards heat flux. Equation (|2ji tracks the burning depth, 
allowing for accretion of new fuel at a rate given by the lo- 
cal accretion rate rii, as well as burning of fuel on a timescale 
E'y./e, where is the energy per gram released in the burn- 
ing. Note that the pressure at the base of the layer v&P = gy 
from hydrostatic balance, where g is the local gravity. We first 
study equations Q and ^ analytically (§2.1), and then show 
some numerical integrations (§2.2). 

2.1. Analytic estimates 

Equations Q and (|2ji constitute a nonlinear oscillator To 
understand its behavior, we consider linear perturbations to 
the steady state solution, which has 



ey = F = rriEi,. 



(3) 



The nuclear energy generation is generally a strong func- 
tion of temperature, and we will assume e ocT", where a = 
dlne/dlnT . The heat flux is approximately F « acT*/3Ky 
(e.g., Bildsten 1998). To simplify the algebra in this section, 
we assume that e depends only on temperature, and that the 
opacity k is a constant. We write the deviations from steady 
state as 6y and ST, giving 
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Defining a thermal timescale for the layer ftherm = CpT /e = 
CpTy/F and an accretion timescale facer = y/ni, and using the 
steady-state relations of equation Q gives 
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I : is useful to write ST = f{t)exp(—t/t; 
§(f)exp(-f /facer)- This simplifies equations ^ and Q, allow 
ing us to combine them into a single differential equation for 
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Equation is the equation for a damped simple harmonic 
oscillator The solution is ST cx exp(A/), with 
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Note that the "damping" term in equation (|8|i can be positive 
or negative depending on the relative temperature sensitivities 
of the heating and cooling. 

The key to understanding the behavior of the burning is 
to note that the timescales ftherm and facer are very different. 
Steady burning at accretion rates near Eddington leads to ig- 
nition at a column depth y w 10^ g cm"^ and temperatures 
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Fig. 1. — Lightcurves for the one-zone model at three different accretion 
rates. The system evolves from the arbitrary initial values v = 2 X 10^ g cm"^ 
and r = 2 X 10* K. At m = 0.95 mEdd. bursts occur with a recurrence time 
of 34 minutes. At m = 0.998 riiEM' oscillations ai'e seen with a period of 
4.7 minutes. At m = 1.05 HiEdd. after a few transient oscillations, the burning 
evolves to a steady state. The steady state flux is mTTnuc, where £nuc ^ 5MeV 
per nucleon. 



r « 5 X 10** K (Schatz et al. 1999). The accretion timescale 
is therefore y/;?? w 1000 s for accretion at m w 10^ g cm"^ s"' 
(roughly the Eddingtonrate). For stable burning F = n'lEi,, giv- 
ing ftherm = cpTy/F = (c pT / E^)(y / m) « 0.01 facer ~ 10 s, since 
for an ideal gas, cpT w {5/2){kBT /nip) = 7.2 keV Tg per nu- 
cleon, whereas « 5 MeV per nucleon (Schatz et al. 1999)"*. 

Because ftherm ^ facer the damping term is usually much 
larger than the oscillatory term, and the system is either pos- 
itively or negatively overdamped depending on the sign of 
(a-4). When a > 4, the damping coefficient is negative, lead- 
ing to strong exponential growth: the steady state solution is 
linearly unstable. When a <4, the damping coefficient is pos- 
itive, leading to strong damping: the steady-state solution is 
linearly stable. Close to marginal stability when a « 4, how- 
ever, the effective thermal timescale ftherm/ |a-4| becomes 
much longer than the accretion time, and the inequality is 
reversed, giving weakly damped or excited oscillations with 
oscillation frequency uj « (2a/faccrftherm)'''^- The condition to 
observe oscillations is that the oscillation frequency should be 

In fact, the gas is partially degenerate at the burning location, with kgT ^ 
Ef, and so there will be a small correction factor to the ideal gas expression 
for £■/>. 
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Fig. 2. — The trajectories in the temperature-coluinn depth plane for the 
models shown in Figure 1 . In each case, the inodel starts at the arbitrary point 
V = 2 X 10** g cm-2 and r = 2 X 10** K. 



larger than the damping rate. 
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Note that this condition can be satisfied when A has either 
positive or negative real part, so that the oscillations may be 
excited or damped. 

Neglecting the slight modification from the damping term, 
the oscillation period is Pose = 2tt/u!, or 
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where we set a = 4 at marginal stability, and take = 5 MeV 
per nucleon. This estimate is remarkably close to the observed 
mHz QPO period of 2 minutes. 

The physics of the oscillations can be understood by con- 
sidering equations (|6} and Q. For most accretion rates, 
the effective thermal time is much smaller than the accretion 
timescale, or ftherm/ |4-a| ^ facer- The perturbations are then 
effectively at constant pressure, or 5y = 0, as commonly as- 
sumed (see, e.g., Fujimoto, Hanawa, & Miyaji 1981; Bildsten 
1998), and equation (|6|l leads directly to exponential growth 
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or decay depending on the relative temperature sensitivities 
of the heating and cooling rates. Close to marginal stability, 
(a-4) Ri 0, however, the effective thermal timescale becomes 
very long compared to the accretion time. Equations ^ and 
O in this limit are 
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These equations nicely summarize the physics of the oscilla- 
tions. Consider an upwards fluctuation in temperature ST > 0. 
At the point of marginal stability, the increase in heating 
rate almost exactly cancels the T'^ dependence of the cool- 
ing rate. The main effect is that the hotter temperature leads 
to faster burning of the accreting fuel, and a decrease in thick- 
ness Sy < on a timescale ~ tg^ca (eq- 1131 ). But a thinner 
layer cools faster since F/y oc l/y^. Therefore, the temper- 
ature fluctuation now begins to decrease, this time on the 
faster timescale ~ ftherm (eq- 1121 ). These changes are out 
phase, driving an oscillation on the intermediate timescale 

~ ('therm'accr) ^ ■ 

The behavior we describe here is shown by the canonical 
example of a nonlinear oscillator, the van der Pol oscillator 
(e.g., Abarbanel, Rabinovich, & Sushchik 1993). This oscil- 
lator consists of an LC circuit with an active element that can 
behave as a "negative resistor", originally a vacuum tube. The 
governing equation is of the foimx+k{x^-l)x+uj^x = 0. De- 
pending on the choice of control parameter k, the behavior of 
this circuit is a limit cycle (relaxation oscillations) with fast 
and slow timescales dominating at different parts of the cycle 
{k> 1), a strongly damped system which evolves to a steady- 
state (A; < -1), or oscillatory with growing or damped oscilla- 
tions (|A:| < 1). These three states are analogous to bursting, 
stable burning, and oscillations near the stability boundary in 
the one-zone model. 

2.2. Numerical integrations 

We have integrated equations ([0 and (|2j in time to deter- 
mine the non-linear evolution of the one-zone model. For the 
nuclear burning, we use the triple alpha reaction (3a — >'^C) 
rate as given by Fushiki & Lamb (1987). We allow for the 
presence of hydrogen in the accreted fuel, however, by en- 
hancing the energy release from the triple alpha reaction by a 
factor Enuc/Eja, where E^a = 0.606 MeV per nucleon is the 
energy release from the triple alpha reaction, and E^uc is the 
energy release from burning the accreted mixture of hydrogen 
and helium to iron group. We assume that Enuc = 1.6 + 4.9Xo 
MeV per nucleon, where Xq is the mass fraction of hydro- 
gen in the accreted layer This expression for finuc includes 
an energy loss of 25% from neutrino emission during ap- 
and rp-process burning (e.g., Fujimoto et al. 1987), and gives 
E,mc = 5 MeV per nucleon for Xq = 0.7, in good agreement 
with the energy release in the steady state burning models of 
Schatz et al. (1999). We write the flux as F = acT'^/3Ky (Bild- 
sten 1998), where the opacity k is calculated as described by 
Schatz et al. (1999). In addition, we include a flux heating the 
layer from below of 0. 15 MeV per nucleon, and a contribution 
to the heating rate from hot CNO hydrogen burning in the ac- 
cumulating fuel layer Neither of these extra contributions to 
the heat balance make a significant difference to our results. 
Note that the amount of hydrogen burned by hot CNO cycle 
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Fig. 3. — The trajectory in the temperature-column depth plane for m = 
0.998mEdd. as shown in Figure 2, but zooming in on the oscillations around 
the steady burning location. 



prior to helium ignition is very small at the rapid accretion 
rates considered here. 

Figure n shows lightcurves from the one-zone integrations 
at different accretion rates. To enable a direct comparison 
with the multizone simulations discussed in §3, we set the 
local gravity to be the Newtonian gravity for a 1.4 Mq, 10 
km neutron star, g = 1.9 x lO'"* cm^ s"'. Throughout this pa- 
per, we define the local Eddington accretion rate to be niEdd = 
8.8 X 10^ g cm~^ s"'. By coincidence'', the stability boundary 
for this one zone model is very close to the Eddington ac- 
cretion rate. In Figure [2 we show lightcurves at ni = 0.95, 
0.998, and 1.05 mEdd- Figure |2] shows the corresponding 
tracks in the temperature-column depth plane. We start the 
simulations with the arbitrary conditions T = 2 x 10*^ K, and 
y = 2 X IQ^ g cm"^. At accretion rates below the boundary, 
the system evolves quickly into a limit cycle corresponding to 
Type I X-ray bursts: slow accumulation of fuel followed by 
rapid burning. Close to the stability boundary, the recurrence 
time is w 30 minutes. At accretion rates above the bound- 
ary the evolution is to a steady state in which the fuel burns at 
rw5x lO'^K and yfulx 1 0^ g cm~^ . For accretion rates very 
close to the transition, m w 1 mEdd, we find oscillations around 
the steady state conditions, with oscillation period w 4 min- 
utes and amplitude a few percent of the Eddington flux. For 
the oscillatory case. Figure |3l shows a more detailed view of 
the trajectory of the solution in the temperature-column depth 
phase space. 

Although the one-zone model is approximate, it is use- 
ful because it allows us to investigate how the properties 
of the oscillations change with parameters such as surface 
gravity and accreted composition. The accreted composi- 
tion could differ from system to system either due to metal- 
licity variations, or variations in the accreted hydrogen frac- 
tion. Intermediate-mass binary evolution models (e.g., Pod- 

^ The physics which stabilizes the burning, the decreasing temperature 
sensitivity of the triple alpha reaction with increasing temperature, has noth- 
ing to do with the physics setting the Eddington luminosity. By coincidence, 
the transition to stability occurs close to the Eddington accretion rate (e.g., 
Bildsten 1998). 
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Fig. 4. — Oscillation period as a function of accretion rate in the one-zone model, for different choices of surface gravity and accreted hydrogen fraction. At 
each accretion rate, we integrate the model for lO'' s, and plot the mean oscillation period, only including data for t > 100 minutes. The solid symbols indicate 
models for which the oscillations grow and reach a steady amplitude; the open symbols indicate models for which the oscillations are damped. 



siadlowski et al. 2002) predict that the companion star in 
many systems is hydrogen deficient, so that the hydrogen 
mass fraction is reduced below the solar composition value 
of Xq k, 0.7. The importance of the hydrogen fraction is 
that the nuclear energy release iSnuc changes significantly with 
only small changes in Xq. In contrast, we expect that metal- 
licity will not have a large effect on the transition accretion 
rate or the oscillation period. This is because the metal- 
licity of the accreted material mainly enters into the one- 
zone model as hot CNO hydrogen burning in the accumu- 
lating layer, but the flux from hot CNO burning, ecNoy ~ 
5 X 10^' erg cm~^ s~' 3'8(Zcno/0.01), is much smaller than the 
steady burning flux, n'lEnuc ~ 5 x lO^"* erg cm~^ s"' (in/mEAd), 
where Zcno is the mass fraction of CNO elements. 

Figure0]shows the dependence of the oscillation period on 
accretion rate for different choices of gravity and accreted hy- 
drogen fraction. We show results for g 14 = 1.9, corresponding 
to the Newtonian gravity of a 1.4 Mq, 10 km neutron star 
(or the general relativistic gravity for a 1.4 M©, 12.3 km neu- 
tron star), gi4 = 2.45, the gravity of a 1.4 Mq, 10 km neutron 
star taking general relativistic corrections into account, and a 
stronger gravity g[4 = 3.1 which corresponds to the general 
relativistic surface gravity for a 2 Mq, 11 km neutron star 
We also show a model with a hydrogen fraction below solar, 
Xq = 0.5. Only the accretion rate range at which oscillations 
are observed is shown. At each accretion rate, we integrate 
the one-zone model for 10* seconds, and plot the mean oscil- 
lation period after discarding the first 100 minutes of data. 

For each choice of ^14 and Xq, the pattern is similar The 
overall range of accretion rates for which oscillations are seen 
is very narrow, Am/m w 1%. For most of this range the oscil- 
lations are decaying. As the stability boundary is approached 
from below, we first see oscillations which reach a steady 



amplitude (indicated by solid squares in Fig. 0} whose fre- 
quency drops rapidly with increasing m. At larger accretion 
rates (open squares in Fig.|4}, the oscillations decay with time, 
on timescales < 1 day, and have a frequency that is less sen- 
sitive to m. The transition from growing to decaying oscil- 
lations is very rapid. In each case, the model indicated by 
the last closed square shows stable oscillations for 10* sec- 
onds, whereas with only a small increment in accretion rate 
the next model (first open square) has oscillations which de- 
cay on a timescale of ^ 10^ seconds. Most interesting is that 
the oscillation period is sensitive to gu and Xq. Increasing 
gravity or decreasing Xq moves the transition from unstable to 
stable burning to higher accretion rates, where the oscillation 
period is shorter. As Xq decreases, the accretion rate range 
over which the oscillations are growing rather than decaying 
is larger 

3. MULTIZONE CALCULATIONS OF BURNING NEAR THE 
STABILITY BOUNDARY 

We now present detailed multi-zone models of nuclear 
burning at accretion rates close to the transition from unstable 
to stable burning. These models are extensions of the calcula- 
tions presented by Woosley et al. (2004) using the implicit ID 
hydrodynamic code KEPLER (Weaver et al. 1978). Woosley 
et al. (2004) calculated sequences of X-ray bursts at accretion 
rates M « 0.03 and 0.1 MEdd- In this paper we show the first 
results of an extension of these calculations to higher accre- 
tion rates. 

Following Woosley et al. (2004), we take the gravitational 
mass of the neutron star to be 1 .4 Mq and R= 10 km, giving 
a Newtonian gravity ^14 = 1.9. The effects of general rela- 
tivity are not included in the simulations itself, but because 
the burning layer is very thin the effects of general relativ- 
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Fig . 5 . — Light curves for different accretion rates (the upper right comer 
gives accretion rate in units of Eddington accretion rate). Panel A shows 
regular bursting with stable recun'ence times. Panel B shows weaker bursts 
with a partial oscillatory behavior in the tail of the burst light curves. Panel C 
shows oscillatory rather behavior and no bursts. Panel D shows very small 
oscillations, essentially stable behavior. The small dips seen in all lightcurves 
on about 15 s time-scale are an artifact of our treatment of accretion, in which 
an entire small outer zone is added periodically. The beginning of this figure 
is 10,000 s after the beginning of the simulation to allow the model to reach 
quasi-stationary conditions. 



ity are small over the extent of the simulated burning layer 
and a local Newtonian frame is a good approximation. Gen- 
eral relativity may be accounted for using appropriate redshift 
factors, as discussed in §4.4 of Woosley et al. (2004). The re- 
sults we provide here do not include these redshift corrections; 
the simulated conditions apply for different combinations of 
neutron star radius and mass that give the same surface accel- 
eration in the local frame, using appropriate scaling of surface 
area, accretion rate, and luminosity. 

Our code includes an adaptive nuclear reaction network that 
automatically adjusts to include or remove isotopes as needed 
to follow the details of the nucleosynthesis, out of a reaction 
rate library of about 5,000 nuclei (Rauscher et al. 2002). The 
calculations presented here use up to 1300 different isotopes. 
The reaction rate library includes recent measurements and 
estimates of critical nuclear reaction rates. The effect of un- 
certainties in these rates are discussed in detail in Woosley et 
al. (2004). The energy generation from the reaction network 
is directly and consistently coupled into the implicit hydro- 
dynamic solver The numerical grid adaptively refines and 
derefines the Lagrangian grid to resolve gradients, but in the 
hydrogen-rich layer in effect is essentially at constant mass 
resolution, i.e., linear in column depth. Accretion is mod- 
eled by periodically adding an extra zone at the surface of 
the star (of column depth w 1.6 x 10^ g cm"^) (Woosley et 
al. 2004). We follow both the compositional and thermal pro- 
files of the layer, including radiative and convective transport, 
with a time-dependent mixing length treatment for convec- 
tion, semiconvection, and thermohaline convection. 

Here we present results for an accreted material metallicity 



of 1/20 solar At an accretion rate M = 0.1 MeiH, Woosley et 
al. (2004) (their model zM) found a sequence of regular bursts 
with recurrence times close to 3 hours. Increasing the accre- 
tion rate in our new sequence of models, we find a transition 
from unstable to stable burning at M = 0.924 Mndd- Figure|5] 
shows the behavior close to the transition accretion rate. The 
lightcurves show a progression from regular periodic burst- 
ing at M = 0.7 Mndd (with recurrence times close to 20 min- 
utes), to a combination of irregular bursts intermixed with 
oscillations at M = 0.923 Mndd, to a regular sequence of os- 
cillations at M = 0.925 Mndd, and finally to stable burning at 
M = 0.95 Medd, with the oscillation amplitude rapidly decreas- 
ing as the accretion rate is increased. 

The oscillation period at M = 0.923 Mem is 1 85 ± 5 seconds. 
Figure |6l shows a portion of the lightcurve at this accretion 
rate. The oscillations have an asymmetric profile, with the de- 
cay lasting twice as long as the rise. Revnitsev et al. (2001) 
noted marginal evidence that the peaks of the mHz QPOs 
were asymmetric, with a steep rise and shallower decline. 
A more detailed comparison of our models with observed 
lightcurves would be interesting, but clearly there should be 
significant harmonic components. The peak-to-peak ampli- 
tude of the oscillation in Figure |6l is « 5 x lO^*^ erg s"', with 
minimum luminosity w 5 x lO"'^ erg s"' and maximum lumi- 
nosity w 10^^ erg s"' . This is in good agreement with the one- 
zone model (compare the middle panel of Fig.^ for a 10 km 
star, a flux of w 10^"* erg cm"^ s"' corresponds to a luminosity 
of sa 10^^ ergs"'). 

The nuclear burning in the oscillation mode is powered 
by ap- and rp-process burning (Wallace & Woosley 1981; 
Schatz et al. 1999), beginning with seed nuclei produced by 
breakout reactions from the CNO cycle, and terminating at 
mass number w 80 (the most abundant nucleus is **"Sr). Fig- 
ure |6l shows the energy generation as a function of column 
depth, radius, and time through several oscillation cycles. Fig- 
ure shows the composition profile of the layer at different 
phases of the oscillation cycle. Beneath the hydrogen burning 
layer, the composition of the ashes shows periodic variations 
with a spacing in column depth of mPosc, where Pose = 185 sec- 
onds is the oscillation period. Note, however, that the hydro- 
gen burning depth is relatively constant during the oscillation 
cycle. The hydrogen burns at a depth « 10^ g cm~^ which is 
« 7 times larger than the amount of mass accumulated in one 
oscillation cycle. The underlying ashes record the variations 
in burning temperature and the resulting variation of the rp- 
process ashes during the oscillation cycle. This is reminiscent 
of the growth of annual rings in a tree trunk. Figure |8l shows 
the distribution of nuclei in the ashes by mass number, and the 
slight variation in the composition through the cycle. 

At the peak in the oscillation light curve, the increased 
temperature in the burning region drives heat transport in- 
wards, heating the underlying material. This is similar to 
the substrate heating during Type I X-ray bursts discussed by 
Woosley et al. (2004), and leads to increase of burning of the 
"^He remaining in the ashes layer by the triple-a process and 
by a captures. In Figure|6]this can be seen as additional spikes 
in nuclear energy generation below the main burning band of 
the hydrogen-rich zone. Of particular interest is the amount 
of carbon remaining in the ashes. Stable burning of accreted 
H/He has been suggested as the source of carbon fuel for su- 
perbursts (in 't Zand et al. 2003; Schatz et al. 2003). Fig- 
ure |S] shows that the amount of carbon at y w 2 x 10^ g cm"- 
is « 2 % by mass, which is very close to the asymptotic value 
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Fig. 5. — Detailed light curve (upper panel) and specific nuclear energy 
generation as a function of time and column depth (lowei panel). In the lower 
thee panels, each darker shading of blue corresponds to a value of energy 
generation one order of magnitude higher; see scale on right hand side of 
the figures. In the second figure from the top we label each depth with a 
Lagrangian column depth. Following a given column depth to the right shows 
the evolution of that fluid element in time. The sloping black line indicates 
the surface of the star (the slope gives the accretion rate). The lower two 
panels show the evolution s a function of radius coordinate. Zero is chosen 
to correspond to the location where hydrogen is depleted. The upper panel 
gives specific nuclear energy generation rate (same as the panel above), the 
bottom panel gives energy generation rate per unit depth. 



of slightly below 2 % in the deeper layers where all the helium 
has been burnt. 

The layering of different compositions in the ashes does 
not persist to great depths. The different layers have differ- 
ent values of the number of electrons per baryon Y^, which 
determines the specific weight of fluid elements under the de- 
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Fig. 7. — Snapshots of structure (temperature: thick gray line; density: 
thick gray dashed line; specific nuclear energy generation: black line) and 
composition (select isotopes, colored lines) during one oscillation cycle; each 
panel is advanced in time by P/4 where P is the oscillation period; the bottom 
panel is advanced by one full cycle. The bottom axis for each figure gives 
column depth, the top axis the corresponding time since accretion began. 
This is the same model as shown in Fig.|5] Panel C. The white and gray stripes 
correspond to one cycle of oscillation each, with the interfaces corresponding 
to the time of a maximum in the light curve at the time of the accretion of that 
layer. The small inserts at upper right corners indicate the position in the light 
curve (red) cycle of the snapshot (black dot; intentionally aligned with a layer 
interface). Note that the decreases of some of the radioactive isotopes in the 
right-hand side of the figure is due to their radioactive decay. An animation 
of this figure can be found at http : / / xraybur st . org/qpo. 



generate conditions in the ashes layer The variation in com- 
position is stable to the Rayleigh-Taylor instability because of 
the thermal buoyancy. Secular doubly-diffusive instabilities, 
however, cannot be suppressed. We find that the thermoha- 
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Fig. 8. — Vai'iation of the composition (summed up for each mass number 
to be invariant against /3+ decays) in the ashes layer at a column depth close 
to 2 X 10** gcm"^ (at the base of the composition profile shown in Fig. 171. 



line or salt-finger instability grows in the layers which have 
an outwardly-decreasing profile. The subsequent mixing is 
apparent in Figure0as the "erosion" of the peaks in the ^'^Ge 
profile and the valleys of the *'Cu profiles at column depths 
above 1 .5 x 10** g cm . The mixing eventually leads to homog- 
enization of the ashes layer at depths > 3 x 10'^ g cm~^. Our 
calculations followed that process to a point where the typical 
column-depth scale of the homogenized regions was several 
times that of the composition oscillations made by the oscil- 
lations in the burning region. 

Schatz et al. (1999) calculated the nucleosynthesis in 
steady-state burning models at m = 1 riiEdd- They found that 
the rp-process endpoint terminated at A w 70, and the carbon 
mass fraction was « 5%. We find slightly heavier rp-process 
ashes (A w 80), and a smaller carbon mass fraction (w 2%). 
This is consistent with the general anti-correlation found by 
Schatz et al. (2003) between the mass of nuclei produced in 
the rp-process and the carbon yield. In our case, a longer 
rp-process gives more time for helium to burn before the hy- 
drogen runs out, leading to less carbon production following 
hydrogen exhaustion. Our burning temperature is approxi- 
mately 10-20% hotter than the models of Schatz et al. (1999) 
which might explain the more extensive rp-process. This is 
perhaps because of a difference in radiative opacities: Schatz 
et al. ( 1 999) use a fit to the results of Itoh et al. ( 1 99 1 ) for free- 
free opacity, whereas the KEPLER code uses the fit of Iben 
(1975) to the radiative opacities of Cox & Stewart (1970a,b). 
We will investigate this difference further in future calcula- 
tions. 



4. DISCUSSION 

The fact that oscillatory burning is naturally expected at the 
transition from unstable to stable nuclear burning was pointed 
out by Paczynski (1981). We have investigated the properties 
of marginally stable burning in this paper with a simplified 
one-zone model (§2) and with detailed multizone simulations 
(§3) using the KEPLER code, extending the Type I X-ray burst 
calculations of Woosley et al. (2004) to higher accretion rates. 
The period, amplitude, and shape of the oscillations agree 
well in both models. Remarkably, the basic physics of the 
oscillations is the same physics underlying a nonlinear relax- 
ation oscillator such as the van der Pol oscillator (e.g., Abar- 
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Fig. 9. — Oscillation frequency as a function of hydrogen mass fraction, 
Xq, for different surface gravities, corresponding to different radii and masses 
of the underlying neutron star. The data is taken from one-zone models 
and rescaled to the oscillation frequency found in the multi-zone model with 
Xq = 0.7 and ^14 = 1.9. Proper general relativistic redshifts corrections and 
radii have been applied. The curve for ^14 = 0.9 corresponds to M = 1.4Mq 
(gravitational mass) and R = 15.2km, the curve for g[4 = 1.48 corresponds 
to M = 1.4Mq and R = 11.9km, ^14 = 1.9 (used in the multi-zone calcula- 
tion) coiTesponds to M = 1.4Mq and R = 12.3km, the curve for gi4 = 2.45 
is for M = 1.4Mq and R = 10km, and the curve for gi4 = 3.1 corresponds 
to M = 2Mq and S = 1 1 km. The dotted lines indicate the boundaries of the 
observed range of QPO frequencies. 



banel et al. 1993). Usually, the positive or negative damping 
term dominates, giving rise to the familiar X-ray bursts at low 
accretion rates or stable burning at a fixed temperature and 
density at high accretion rates. Close to the marginally stable 
point, however, the effective thermal timescale is very long, 
allowing the underlying oscillation period of the system to be 
seen. This period is close to the geometric mean of the ther- 
mal time and accumulation time of the burning layer (eg. II II . 

This behavior naturally reproduces two properties of the 
mHz QPOs observed by Revnivtsev et al. (2001): the ob- 
served periods of w 2 minutes, and the fact that mHz QPOs 
were observed in only a narrow range of luminosities in 
4U 1608-52, (0.5-1.5) x 10^^ erg s"'. Identification of the 
mHz QPOs with marginally stable nuclear burning would for 
the first time relate a feature of the persistent X-ray emission 
to the neutron star surface in sources which are not X-ray pul- 
sars. Further, our one-zone models indicate that the oscil- 
lation period is very sensitive to the surface gravity and ac- 
creted hydrogen fraction (Fig. |4}. One of the difficulties in 
comparing X-ray burst properties with theoretical models is 
the uncertain relation between X-ray luminosity and accre- 
tion rate (e.g., Cumming 2003). This uncertainty is removed 
for marginally stable burning because it occurs at a specific 
accretion rate. The dependencies on surface gravity and hy- 
drogen fraction need to be confirmed with multizone models, 
although unfortunately scanning the parameter space for the 
location of the transition is computationally very expensive. 
Our one-zone model results (Fig.|4} suggest that the observed 
periods of « 2 minutes require either < 0.7, as predicted 
by intermediate mass evolution models (e.g., Podsiadlowski 



9 



et al. 2002), or surface gravity gi4 w 3, for example corre- 
sponding to a 2 Mq (gravitational mass) star with R= II km. 
This can be seen in Figure |9] which shows the oscillation 
period as a function of Xq and gi4. In this figure, we have 
rescaled the one-zone model results to match the multizone 
model for Xq = 0.7 and g[4 = 1.9; in addition, we include the 
gravitational redshift factor^. Future comparisons of theoret- 
ical models with mHz QPO periods and lightcurves are po- 
tentially sensitive probes of the surface gravity and accreted 
composition. 

Two observed features of mHz QPOs remain to be ex- 
plained. The first is the Q value of the oscillation, which 
Revnivtsev et al. 2001 found to be g = v j « 3-4. This 
may be related to the range of accretion rates for which os- 
cillatory nuclear burning can be observed. Revnivtsev et 
al. (2001) constrained the range of luminosities at which mHz 
QPOs are present to be (0.5-1.5) x lO" erg s"' in 4U 1608- 
52. In contrast, the theoretical range where oscillations are 
seen is much smaller, within a range /S.in/m k, 0.01 around the 
transition accretion rate. Moreover, for much of this range, the 
oscillations decay in time. Further observations which con- 
strain the range of luminosities for which mHz QPOs can be 
observed would be valuable. 

The second puzzle is that our theoretical models, in agree- 
ment with previous estimates (Fujimoto, Hanawa, & Miyaji 
1981; Ayasli & Joss 1982; Bildsten 1998), find that the tran- 
sition to stable burning occurs at an accretion rate close to the 
Eddington rate (M = 0.924 M^id for the model presented in 
§3). In contrast, the X-ray luminosity at which the mHz QPOs 
are observed in 4U 1608-52 is w 0.5-1.5 x lO" erg s"', im- 
plying an accretion rate ten times lower, M w 0.1 MEdd- The 
same factor of ten appears when comparing the theoretical and 
observed oscillation amplitudes. The observed amplitudes of 
mHz QPOs correspond to flux variations of « 1-2% (Revnivt- 
sev et al. 2001). As noted by Revnivtsev et al. (2001), this is 
consistent with the ratio of nuclear energy release from burn- 
ing the accreted material to gravitational energy release from 
accretion. We find fractional amplitudes of the expected few 
percent level in the theoretical models. The absolute peak- 
to-peak flux variation is a factor of ten larger than observed, 
however, roughly 5 x lO^-' erg cm"- s"' or 5 x lO-'^ erg s"'. 
This is more than 30% of the observed persistent luminos- 
ity. Note that the oscillation amplitude decreases rapidly with 
increasing accretion rate above the transition accretion rate, 
so there does exist a small range of accretion rates where the 
amplitude does match the observed value. 

A simple explanation for both of these discrepancies is that 
the local accretion rate onto the star is close to the Edding- 
ton rate, even though the global accretion rate is much lower. 
Because the burning layer is very thin, the properties of the 
burning depend only on the local accretion rate, which may 
vary across the surface of the star. In particular, if the ac- 
creted material spread over only sa 10% of the surface, the 
local accretion rate would be ten times larger, comparable to 
the Eddington rate, and the emitted luminosity would be ten 
times lower than if the burning covered the entire stellar sur- 
face. This would allow marginally stable burning to occur at 
a global rate of « 0.1 Eddington, while giving flux variations 

^ Gravitational redshift increases the oscillation periods shown in FigureRI 
by 30%. The numerical results in §3, however, suggest that the one-zone 
model overpredicts the oscillation period by a similar factor. Therefore, the 
oscillation periods from the one-zone model without redshifting are in fact 
approximately those we expect from multizone models corrected for gravita- 
tional redshift. 



of a few percent as observed. 

The physics that might cause confinement of the fuel onto 
a small fraction of the surface is not obvious. The pressure at 
the base of the burning layer \s P = gy= 10^^ erg cm"-' giAy%, 
implying that magnetic fields of strength approaching w 
VSttP w 3 X 10" G would be required to confine the fuel. 
This is much larger than the « 10^-10^ G fields assumed for 
the neutron stars in LMXBs (believed to be the progenitors 
of the millisecond radio pulsars; Bhattacharya et al. 1995), 
although small scale fields of these strengths might exist on 
the surface. The need to transport angular momentum could 
also potentially delay spreading of material accreted from a 
disk onto the equator of the star. Inogamov & Suny aev ( 1 999) 
studied this problem with a one-zone model of the spreading 
layer and a basic prescription for angular momentum trans- 
port. They found column depths < 10"* g cm"- in the spread- 
ing layer, much smaller than the burning depth. 

The possibility that the covering fraction changes with ac- 
cretion rate was suggested previously by Bildsten (2000), 
but in the opposite sense, with covering fraction increasing 
with m. The motivation was to explain a puzzling change 
in burst behavior that is observed to occur at a luminosity ~ 
10^^ erg s"'. EXOSAT observations of several Atoll sources 
showed that as X-ray luminosity increased, burst properties 
changed from regular, frequent bursts (frecur ~ hours) with 
energetics consistent with burning all of the accreted fuel in 
bursts, to irregular, infrequent (fi-ecm- > 1 day) bursts whose en- 
ergetics indicate that only a small fraction of the fuel burns in 
bursts (van Paradijs et al. 1988). RXTE and BeppoSAX ob- 
servations confirmed this result, with particularly good cover- 
age for the transient source KS 1731-260 (Muno et al. 2000; 
Cornelisse et al. 2003). Cornelisse et al. (2003) found that 
observations of nine bursters with BeppoSAX were consis- 
tent with this pattern of bursting, with a universal transition 
luminosity Lx ~ 2 x lO^*^ erg s"'. 

This change in bursting behavior is not predicted by the 
standard theory, in which regular bursting should continue 
up to the stability boundary at n'l « WEdd- Several theoreti- 
cal explanations for the discrepancy were put forward, includ- 
ing mixing by Rayleigh-Taylor (Wallace & Woosley 1984) or 
shear instabilities (Fujimoto et al. 1987) which might allow 
more rapid burning of hydrogen, a new mode of burning in- 
volving slowly propagating fires at n'l > 0.1 WEdd (Bildsten 
1993), or that the covering fraction of accreted material in- 
creases at higher accretion rates, lowering the accretion rate 
per unit area and lengthening the time between bursts, giv- 
ing hydrogen time to burn stably (Bildsten 2000). We also 
mention here that Narayan & Heyl (2003) calculated linear 
eigenmodes of truncated steady-state burning models, and 
found stability for accretion rates M > 0.25 MEdd, more con- 
sistent with observations. The lower accretion rate is likely 
the cause of the small oscillation frequencies that they found 
for marginally stable burning (period ~ 20 minutes). We find, 
however, that bursting continues unabated up to M « MEdd- 
Further work comparing linear stability analysis with numer- 
ical calculations seems to be required. 

The observations of mHz QPOs at a luminosity close to 
w 10^^ erg s"' and their interpretation as marginally stable os- 
cillatory burning at a local accretion rate rii w mEdd provide 
new input for these ideas. As we have discussed, if changes 
in the covering fraction are responsible, this suggests that the 
covering fraction decreases rather than increases at this lumi- 
nosity. The observed Type I bursts at Lx > 10^^ erg s"' could 
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be accommodated if there was a slow "leak" of fuel away from 
the stably burning region. This fuel would deplete hydrogen 
as it accumulated, giving occasional short helium-rich bursts. 
Other mechanisms, such as stable burning driven by mixing of 
fuel by shear instabilities might also lead to oscillatory burn- 
ing. More theoretical work is needed. One clue is that the 
ultracompact source 4U 1820-30 which most likely accretes 
pure helium (Bildsten 1995; Gumming 2003) shows a similar 
transition at a similar luminosity, implying that the nature of 
the transition does not depend on accreted composition. 

This question is also likely to be relevant for superbursts 
and Type 1 burst oscillations. Superbursts are long duration, 
rare, and extremely energetic Type 1 X-ray bursts (up to 1000 
times the duration and energy and less frequent than normal 
bursts) believed to be due to unstable carbon ignition (Gum- 
ming & Bildsten 2001; Strohmayer & Brown 2002). Super- 
bursts are only observed at luminosities above w lO^*^ erg s~', 
and from sources for which burst energetics indicate that 
bursts burn only a small fraction of the accreted fuel (in 't 
Zand et al. 2003). This fits nicely with the theoretical result 
that stable burning is much more efficient than unstable burn- 
ing at producing the carbon fuel (Schatz et al. 2003). How to 
achieve stable burning theoretically at m < liiEdd, however, has 
been an open question. Type 1 burst oscillations are high fre- 
quency oscillations during Type 1 X-ray bursts believed to be 
due to burning asymmetries on the surface. Burst oscillations 
are preferentially seen at higher accretion rates, in the banana 
branch of the color-color diagram (e.g., Muno et al. 2000). 
Incomplete covering of the surface might help to explain the 



origin of burst oscillations, facilitating inhomogeneous burn- 
ing when Type 1 X-ray bursts are able to occur 

Finally, we note that Atoll sources undergo a transition from 
the island state to banana branch in the color-color diagram at 
a luminosity of Lx ~ 10^^^ erg s"'. It is well-known for these 
sources that X-ray luminosity does not track accretion rate on 
short timescales (van der Klis 2001). One explanation for the 
transition from island state to banana branch in Atolls is that 
a hot quasi-spherical accretion flow at low rates is replaced 
by a thin disk at high rates (e.g., Gierlinski & Done 2002). If 
this picture is correct, it could be that the change in accretion 
geometry affects the distribution of fuel on the neutron star 
surface. An interesting observational question is whether the 
appearance of the mHz QPOs is linked to a particular lumi- 
nosity range, or a particular part of the color-color diagram, 
for example the island to banana transition. 
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